#include "slepdef.h"
#include "slepblas.h"
#include <string.h>
#include <math.h>
#include <assert.h>


void slep_dgemm(slep_transposity transA, slep_transposity transB, double alpha,
				slep_matrix* A, slep_matrix* B, double beta, slep_matrix* C){
	int i,j,k;
	double temp;

	// For now we implement only for dense matrices
	assert( A->sparsity == SlepDense && B->sparsity == SlepDense && C->sparsity == SlepDense);
	assert( (transA==SlepTranspose || transA==SlepNoTranspose) && (transB==SlepTranspose || transB==SlepNoTranspose) );
	if(transA == SlepNoTranspose && transB == SlepNoTranspose ){
		assert( (A->n==B->m) &&(A->m==C->m)&&(B->n==C->n));
		for(i=0; i< A->m; ++i)
			for(j=0; j< B->n; ++j){
				C->data.mtrx[j*C->ld + i] *= beta;
				temp = 0.0;
				for(k=0; k < A->n; k++) // A->n == B->m
					temp += A->data.mtrx[k*A->ld + i] * B->data.mtrx[j*B->ld + k];
				C->data.mtrx[j*C->ld + i] += alpha * temp;
			}
	}else if(transA == SlepNoTranspose && transB == SlepTranspose ){
		assert( (A->n==B->n) &&(A->m==C->m)&&(B->m==C->n));
		for(i=0; i< A->m; ++i)
			for(j=0; j<B->m; ++j){
				C->data.mtrx[j*C->ld + i] *= beta;
				temp = 0.0;
				for(k=0; k < A->n; k++) // A->n == B->n
					temp += A->data.mtrx[k*A->ld + i] * B->data.mtrx[k*B->ld + j];
				C->data.mtrx[j*C->ld + i] += alpha * temp;
			}
	}else if(transA == SlepTranspose && transB == SlepNoTranspose ){
		assert( (A->m==B->m) &&(A->n==C->m)&&(B->n==C->n));
		for(i=0; i< A->n; ++i)
			for(j=0; j<B->n; ++j){
				C->data.mtrx[j*C->ld + i] *= beta;
				temp = 0.0;
				for(k=0; k < A->m; k++) // A->m == B->m
					temp += A->data.mtrx[i*A->ld + k] * B->data.mtrx[j*B->ld + k];
				C->data.mtrx[j*C->ld + i] += alpha * temp;
			}
	}else if(transA == SlepTranspose && transB == SlepTranspose ){
		assert( (A->m==B->n) &&(A->n==C->m)&&(B->n==C->n));
		for(i=0; i< A->n; ++i)
			for(j=0; j<B->m; ++j){
				C->data.mtrx[j*C->ld + i] *= beta;
				temp = 0.0;
				for(k=0; k < A->m; k++) // A->m == B->n
					temp += A->data.mtrx[i*A->ld + k] * B->data.mtrx[k*B->ld + j];
				C->data.mtrx[j*C->ld + i] += alpha * temp;
			}
	}
}


void slep_dsymm(double alpha, slep_matrix* A, 
				slep_matrix* B, double beta, slep_matrix* C){
}


void slep_dtrmm(slep_side side, slep_transposity transA, double alpha, 
				slep_matrix* A, slep_matrix* B){

}



void slep_dtrsm(slep_side side, slep_transposity transA, double alpha, 
				slep_matrix* A, slep_matrix* B){

}